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I. Introduction- 

The problem of accurately modeling the relative positions of two orbit- 
ing vehicles has become increasingly important with the advent of missions 
whose success depend upon this information.. Three instances of* such missions 
are: 

1) Missions involving rendezvous 

2) Missions which attempt to study certain physical phenomena by observ- 
ing their relative effects on two orbiting vehicles. 

3) Inexpensive modeling of a set of satellites to check for possible 
col 1 i si on . 

One obvious approach to determining the relative position of two vehi- 
cles is simply to difference the position vectors of the two vehicles. Thus, 
if the position vectors, F.j and of the two vehicles are given as: 

A A /V 

r.| = i + y^ j + Z.J k (1) 

and' 

^2 ~ ^2 ^2. ^ ^2 ^ * 
the relative vector, r, would be computed as 

F = Fj “ ^2 ^^1 " ^ 2 ^’ ‘‘‘ ^^1 ’ ^ 2 ^ ^ ^^1 ■ 

Clearly, when r-^ approaches Fg, a serious loss of significant digits 
may result from the subtractions involved due to the finite word length 
of the computer being used. 

There are .obviously other limiting factors upon the accuracy with which 
the relative vector can be obtained. The most visible limiting factor if the 
method used to obtain the' positions of each- satellite, in the-case of -numeri cal 
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integration this is a factor of the accuracy of the integrator, the stepsize 
used and the length of 'time integrated. 

The use of time transformations, have recently been applied to the equa- 
tions governing the motion of space vehicles in attempts to more accurately 
predict their positions. It is hoped that these methods may be benificial in 
the case of relative positions between two satellites. These studies have 
shov/n‘ that the value of the exponent used in the time transformations, n, can 
have dramatic effect on the global integration error. These investigations 
have generally been concerned with applying the transformations to a single 
orbiting body. Velez^^ has shown that the transformations with 1 < n < 2 
enhance the stability of the system in the sense of Liapunov, and that as n 
is increased, the numerical stability of the numerical integration algorithms 
increases. In addition, Velez^^’^^ Feagin and Mikkilineni^, and Nacozy^ 
have studied the effects of these methods on truncation errors. Nacozy® gives 
an excellent survey of these investigations. One major result of these . 
studies is that the optimal value of n is dependant upon the model employee 

Other studies such as Dunning^ and Meuller^*^ haye developed special- 
ized equations of motion based upon approximations which are valid only 
when the vehicles are in close proximity. For this reason, when using 
these rendezvous type equations it is necessary to alternate between a 
set of equations which are valid only for near approaches, and a set valid 
for large separations. Such switching can introduce a number of complica- 
tions. One primary complication introduced is the determination of when 
to switch, which increases the complexity of the required coding. Other 
complications resultsfrom questions inherent in the alteration, such as 
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how to restart the integration if a multi-step integrator is used, and how 
to control loss in accuracy resulting from required transfonnations. 

This study compares eight formulations of the equations of motion 
for two bodies orbiting- a central body. All equations compared will be 
valid for both close approaches and large separations. As well as avoiding 
the problems mentioned above, this will allow a more careful study of 
the regions in which the individual equations are most applicable. 



CHAPTER II 


Formulations Studied 

Eight equation sets are studied which represent several philosophies 
for improving computational accuracy in the relative motion problem. Some 
of the formulations attempt to increase precision by coupling the equations. 
This is accompl-i shed by integration of one satellite in earth centered 
coordinates and integration of the relative vector joining the satellites. 
Another concept employed attempts to make the equations of motion more 
stable through the use of a new independent variable of integration in order 
to reduce the propagated error. Several studies have shown that time 
smoothed equations enhance the stability of the equations, [Beaudet * ], 
[Velez^^’^^]. Still other formulations are based upon a more direct 
approach and attempt to reduce the errors resulting from subtraction of the 
two nearly equal quantities which appear while computing relative acceler- 
ation. The formulations studied are listed below. 

1. Standard Cowell 

2. Cowell equations modified for relative motion 

3. Independently time smoothed 

4. Time smoothing of one satellite and the relative vector, based 
upon smoothed satellite 

5. Time smoothing of one satellite and the realtive vector, based up 
upon both satellites 

6. Cowell equations modified for relative motion employing Nacozy- 
Szebehely method 
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7. Nacozy-Szebehely method applied to equation set Four 

8. Nacozy-Szebehely method applied to equation set Five 

The most cotimonly used formulation of the equations of motion in space 
dynamics is the Cowell set. Two versions of these equations were included 
in the comparisons to serve as reference. 

2.1 The Cowell type equations 

Equation Set 1 (Standard Cowell equations,) 

The classical Cowell equations were applied independently to each 
satellite as follows: 

Ti = - u ^ + P (r^, Fj) 

— ~ • (3) 

*■2 ° (>' 2 > 

where the first term on the right is Newtonian Force and the P‘s represent 
the perturbing forces. 

The relative vector r is determined at desired output points by 
differencing 

( 4 ) 

In this formulation the equations are very easy to implement, however, 
the problem of lost significant digits in the computation of r from (4) 
can lead to difficulties in obtaining the relative vector accurately when 
r^ and are very close to one another. In an attempt to reduce this 
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source of error, it is shown in the next section that the equations 
can be coupled by altering them for integration of one satellite re- 
ferenced to the central body, and the other satellite position rela- 
tive to the first; 

Equation Set 2 . (Cowell .Equations modified for relative motion) 
Differentiation of equation (4), and substitution of equations (3) 
into the result yields 



In this formulation, one of equations (3) and equation (5) are integrated 
simultaneously. In the evaluation of equation (5)', equation (4) must be 
used to determine the position vector of the satellite which is not integrated 
in earth centered coordinates. It should be noted that as approaches r^, the 

subtraction required by equation (5) may still result in a loss of signifi- 
cant digits. Thus, the original problem has not been eliminated. Rather 
than having the loss of significant digits appear directly in the calcula- 
tion of F, as F = F^ - ?2, it now appears in the acceleration of F. 

There are other techniques which can be applied to the problem which may 
allow the solution to be computed with greater precision. One such tech- 
nique is the use of time smoothing to enhance the stability of the equation. 

2.2 The Time Smoothed Equations 

The time smoothed eaua+ions are all based on a transformation of the 
indepe.ndent variable of integration from time to a new independent variable 
s. This transformation generally requires that another differential 
equation, relating t and s, be added to the equation set. In 
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general, the-'irelation can be given as 
ds = f dt, 

where f is a function of the magnitudes of the position vectors r^ 
and/or Fg. Detailed derivations of the time smoothed equations are 
presented in the Appendix. The most obvious use of these equations is to 
apply them independently to each satellite. 

Equation Set 3 ; (Independently smoothed equations for each satellite) 
For this set f^ and fg are defined as • 


and 


ds., 1 

^ dt 

ds2 1 


The resulting equations become 


(7) 

( 8 ) 


dr^ dr 


1 


ds-| ds.| 


T ^2 y. 2n,-3 - . ^ \ ^2 2n-, 


(9) 




where primes denote differentiation with respect to 

variable of integration. The relations between s^, 

^It n, 
a^ r^ ds^ 



s 


10 



ds 


2 


s, the new independent 
$ 2 * and t are 


(10) 


and 
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Thus, in order to obtain ihfonnation about the relative vector at a 
desired output point it is first necessary to match t^ and tg. This time 
matching is extremely difficult since the relations are both highly non- 
linear, and get even more so as the values of n^ ana ng increase. For 
this reason a formulation which eliminates the need for this time matching 
is desirable. One method for accomplishing this is integration of the 
relative vector directly. 

Equation Set 4 : CTime smoothing based on one satellite's distance from the 
central body applied to relative vector.) 

One satellite and the relative vector are smoothed, the smoothing 
function f is defined as 


4 T _ ds 
^ ■ dt 


1 




(n) 


The resulting equations of motion are: 


Hi HliHli . „ .2 r + F ? 1.2 r2n? 

ro ds ds ^ ^2 ^2 ^ r,j 2 


2 ” 2'“2 2 


_ ^2 dr dr 
rg ds ds 


.2 ,.2n2 


+ ag rg P (r^,r^)-P(rg,rg)+{i ^ 2 


In order to evaluate (13) the relation 


02 ) 

(13) 


r^ = rg + r 


(14) 


must be used. For this formulation, there is only one relation between t 
and s, . 




(15) 
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This eliminates the major drawback inherent in equation set 3, the required 
matching of- t^ and t 2 .’ However, only information about one of the 
satellites is being incorporated into the smoothing function. If the 
smoothing function is based upon r-j and rg, it might adapt the equations 
to the satellites even more accurately. 

Equation Set 5 : (Time smoothing based on both satellites applied to the 
relative vector) 

Another time smoothing formulation, suggested by Dr. V. 6. Szebehely®, 
also smooths one satellite and the relative vector, but this time 
information about both satellites is supplied to the smoothing function 
which is now defined as 




~n . 


(16) 


For this case the equations of motion emerge as: 


_ H n _ 2n rp _ _ ^ 2n 

'-2 = KTT -3 * ('-2>'"2>(''i'’2) 

1 Z rp 


07) 


r = 




(r,r2'+r2r,')r'-u(rir2> 


2n 


ri r2 

3 " 3 

+ 

• • 

ri \ 




<'-r2> 


2n 


The relationship between t and s is 



( 18 ) 



10 


In each of the time smoothed equation sets , the exponent in-.' the time 
smoothing equation, denoted n, is arbitrary. The value of n can 
affect the methods to a great extent, as will be shown. Thus, one 
necessary consideration of the present study is to determine the most 
advantageous choice for n for each time smoothed equation set. In this 
study, each equation set was tested at n=1.0,- 1.5, and 2.0. 

It can be seen that in each of the formulations involving integration 
of the relative vector, there is a term arising in the acceleration of r 
due to the two body forces of the form 



As discussed earlier, this subtraction can lead to numerical difficulties. 

The Nacozy-Szebehely method described below attempts to relieve these 
difficulties by computing this term in a more suitable fashion. 

2.3 The Nacozy-Szebehely method for computing relative acceleration. 

O 

Na cozy- and Szebehely have recently developed methods for accurately 
calculating the acceleration of the relative vector. Nacozyand Szebehely 
apply Encke's method as modified by Potter to the problem of non-interacting - 
space vehicles. 

Expressing the relative vector, r as 



(19) 
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the relative acceleration can be written in the following form. 


r = V 



+ Hryr^) 


_ 

“ P 


( 20 ) 


The first tenn on the right hand side is the acceleration due to the 

__ • __ __ _• 

Keplerian forces, while-. P(r.|,r.j) and P {^ 2 *^ 2 ^ perturbing forces. 

If the perturbing forces are small, compared to the Keplerian forces, the 
errors incurred in calculating the Keplerian acceleration can be expected 
to dominate. In cases where the perturbing forces are large, a similar 
technique for calculating their differences as described below for the two 
body forces may be employed. 

Observing that in equations (20), the Keplerian term is similar to a 
term arising in Encke's method for special perturbations, Nacozy and 
Szebehely applied a technique analogous to Encke's to that term. Rewriting 
the subtraction in equation (20) as 



a quantity q is defined as 


(rg + V) + r 


( 21 ) 



as in tne classical Encke method. It should be noted that when q is small 
(the magnitude-of F.j is close to the magnitude of r2) , a substraction of 
two nearly equal numbers is still involved. A more accurate means of 



computing q can be obtained by substitution of equation 09 \ into the 
above expression. This results in 

~( 2 r 2 +r) * r 

(Fg+r) (r 2 +V) (22) 

The classical approach at this point is to define 



3/2 

(1 + q) - 1 


(23)' 


and expand f(q) in a power series in q,- or to use tabulated values. 

Batting however, derives a closed form expression for f{.q) that pro- 
vides accurate computations of f(q) for small q. 

Nacozy and Szebehely adopt this formulation, which is designed for 
mutually attracting bodies to the relative motion on non-attracting bodies. 
Multiplying equation (23) by 

1 + (1 + 

1 + (1 + q ) 2/2 


equation (23) becomes 
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Substitution into equation (19) for the relative acceleration yields; 



+ (1 + f(q))r + P(r 


]■ 


1 


.ri) 


- P(r2,r2^ 


(25) 


This method for computing the relative acceleration was incorporated 
in the Cowell and Time smoothed equations involving integration of the 
relative vector. 

Equation Set 6; (Nacozy-Szebehely (N-S) method applied to the modified Cowell 
Formulation) 

For this formulation, the- equations for Set Z are modified to use the 
(N-S) technique for computing the relative acceleration. The equations 
become 

r 

+ p (F2, ■ 

r ^ [f(q) + (1 + f(q))^^ + P(F,,r^) -P‘(F2,72) 

The implementation of the N-S method into the rime smoothed equation sets 
is straight forward. The N-S formulation for the acceleration of the relative 
vector can be substituted directly where it is required. 

Equation Set 7 ; (Nacozy-Szebehely method applied to Equation set 4). 

When the N-S method is substituted into equation set 4, (equations (12) and 
(13)), the results are. 
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II 


_ji 

r 



+ 1 + 


dro dr. 


2 

ds 


'ST “ ^ ® 2’^2 


2^ 2n-3 


^2 + P(Kj r^)“ P 





(27) 


The Nacozy-Szebehely method can be incorporated into equation set 5 
in the same manner: 

Equation Set 8: (N-S method applied to equation set 5.) 

After introduction of the N-S technique into equation set 5, the equations can 
be written as 


''2" = "F^ ‘'■r2'-"''2''l')'-2' - M P 

n , 2n _ _ ^ ^ . (28) 

+ + 0 + m)7]j 


Eight equation sets are now available f or .compari son , the two Cowell 
formulations, the three time smoothed equation sets, and the Nacozy-Szebehely 
method applied to Modified Cowell formal tfl ati on , and the two time smoothed 
eouations involving integration of the relative vector. 

The next consideration is how the comparison was carried out, the 
first aspect being which satellite pairs were modeled. This is described 
in the next chapter on Computational Procedures. 



CHAPTER III 


Computational Procedures 

Obviously the eight equation sets cannot be compared for every con- 
ceivable pair of satellites. Thus, a representative group was chosen which 
includes many different situations. The effect which varying the altitudes 
of one or both satellites has on the relative efficencies-of the equation 
sets oan be studied as well as- the effects of the eccentricities of the 
orbits. 

Due to the force model chosen (which will be described In the next 
section of this chapter), the absolute inclination, (I) ascending node (to), 
and argument of perigee («) have no effect on the models. Only the relative 
inclination (i), eccentricity (e), and semi -major axis (a) of the two 
orbits will affect them. In all cases inclination to and were 
chosen as 0. In all pairs, closest approach is at perigee and the initial 
time 2,000 seconds after the start of the integration. This allows the 
first small separation to have an effect on the integration. In an attempt 
to observe the effects of the altitudes of the orbit on the equation sets, 
very high, moderately high and near earth orbits were modeled. The effect 
of the subtraction on the equation sets was studied both by using orbits 
which were similar and dissimilar as well as by using orbits which were 
close for the entire integration time and those which were close only for 
small periods of time. 

3.1 Satellite Pairs Modeled 

The eight satellite pairs modeled are diverse enough to cover a large 
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variety of possible orbits while they have enough similarities to allow 
comparisons of the effects of changes in certain parameters on the relative 
efficiencies of the methods. 


Pair A. 

This first pair was inspired by the International Sun-Earth Explorers 
(ISEE) A B satellites which will be launched to study the structure of the 
solar wind. This pair models two high altitude orbits which have large 
eccentricities. The orbital parameters are: 


Semi -major axis 

Eccentricity 

Inclination 

Peri focus passage 

Argument of perigee 

Ascending node 

Final integration time 


Satellite 1 
76,000km 
.913 

QO. 

2,000 sec 
0 ° 

0 ® 

150,000 sec 


Satellite 2 
76 ,000km 
.78142 
0 ° 

2,000 sec 
0 ® 


Figure la shows these orbits, while Figure lb displays the time history 
of the relative distance between the two satellites. The minimum separation 
for this pair 6ccur at apogee and perigee and are on the order of 
10,000 km. Apogee for satellite one is 145 380km, and for the second 
satellite is 135 380 km, so for this pairs the subtraction of 
at its worst case involves only minimal loss in significance. This allows 
a comparison of the equation sets when the subtraction does not introduce 
a significant error. 

The next pair modeled is more representative of the ISEE A B pair, with 
much smaller separations at apogee and perigee to determine the relative 
effects on the equation sets of the close approach. 
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Pair B 


The orbital parameters for this pair are: 

Satellite 1 


Semi -major axis 
Eccentricity 
Inclination- 
Peri focus passage 
Argument of perigee 
Ascending node 

Final integration time 


Satellite 2 
76,000km 
.9129342105 
0 ® 0 ® 

2,000 sec 2,000 sec 

0 ° 0 ° 

0 ® 0 ° 

150,000 sec. 


76,000km 

.913 


These parameters describe orbits which have minimum separations of 
only 5 kilometers, thus at agogee and perigee, as many as five significant 
digits can be lost in the subtraction of ^2 * 

Figure 2a shows the similarity between the two orbits (one cannot 


distinguish between them at the scale used) while 2b displays the relative 


separation as a function of time. The close approach is the primary 
difference between this and the previous pair, and will allow us to draw 
conclusions as to the effect of a very small separation for this type of 
orbit. 


The effects of these large eccentricities on the equation sets is 


also of interest. One way to study this is to use another pair model with 
much smaller eccentricities. 
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Pair C 

These satellites are two high altitude orbits, similar to the previous 
set, but one is circular, while the other is very nearly circular. 


Orbital Parameters 
Satellite 1 

Semi -major axis 42,162.0km 


Satellite 2 
42,162.8km 


Eccentri ci ty 
Inclination 


0 l.OE-6 

0 ® 0 ° 


Peri focus Passage 


2,000 sec 


2,000 sec 


Argument of Perigee 0® 

Ascending Node 0® 


Final Intigration time 


86,160 sec 


0 ° 

0 ® 


As can be seen from Figure 3a, these orbits are very similar, as were Pair 

-2 

B. The minimum separations (from Figure 3b) are approximately 4.216x10 km., 
and again up the five digits can be lost in the subtraction of Fg “ 

Thus, the major difference between Pairs B and C are the eccentricities of 
the orbit. 

Up to this point, only very similar orbit pairs have been described 
but the comparisons of the equation sets for dissimilar orbits should 
also be investigated. 


Pair D 

For this case two very dissimilar orbits are modeled, one circular 
and the other moderately eccentric with close approach at perigee. The 
period of the second satellite is three times that of the first, so that 
every three revolutions of the first satellite there is another close 
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Jt vs. y for Pair C 
(kiloreters) 



0.00 20.00 40.00 60.00 80 00 100.1 

TIME *10"^ seconds 


Figure 3b 

TIME HISTORY OF RELATIVE DISTANCE FOR 
PAIR C 


OEIGINAL PAGE IS 

OP POOR quality: 


22 



X V5. y for Pair » 
(kllonwtcrt) 


ORIGMAi; PAGE IS 
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approach. The orbits and their relative time histories are shown in 
Figures 4a and 4b. 
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Orbital 

Parameters 

Satellite 1 

Satellite 2 

Semi -major axis 

42.162.8 km 

66,929.37106 

Eccentricity 

0 

.37004010125 

Inclination 

0® 

O'* 

Peri focus passage 

2,000 sec 

2,000 sec 

Argument of perigee 

O'* 

O'* 

Ascending node 

0° 

0° 

Final integration time = 

86,160 sec 



At the point of minimum separation, up to five significant digits may 
be lost in computing r^ " ^2* close separation occurs only once 

every three revolutions of Satellite 1, we can attempt to observe the 
propagated effects of the subtraction of the two nearly equal numbers. 
Another objective of the study is to compare the relative efficiencies 
of the equation sets as the altitudes of the orbits vary. 

Pair £ 

This pair models two near earth low eccentricity orbits which are 
very similar and have small separation for the entire orbit. The only 
difference between this and Pair C are the altitudes of the orbits 

The orbits are indistinguishable at the scale shown as can be seen in 
Figure 5a. The relative distance time history is shown in Figure 5b. 
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Orbit Parameters 


Semi -major axis 

Eccentricity 

Inclination 

Peri focus passage 

Argument of perigee 

Ascending Node 


Satellite 1 
7378.165km 

0 

0 ® 

2,000 sec 
0 *^ 

0 » 


Final Integration Time = 10,000 sec 


Satellite 2 . 
7378.165 km 

1 .OE-6 

2,000 s>ec 
0 ° 

0 ° 


At the two points of minimum separation approximately 5 digits are 
lost in the subtraction of the two position vectors. 

To further study the effects of eccentricities on the equation sets 
a situation similar to Pair E is used, but this time with a greater eccen- 
tricity for both satellites. 


Pair F 

Pair F represents two low altitude orbits which are very similar, 
but of moderate eccentricity. The orbit shape and relative vector time 
history are shown in Figure 6a and 6b, respectively.. 
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Orbit Parameters 


Semi -major axis 

Satellite 1 
- 7378.165km 

Satellite 2 
7378.165km 

Eccentricity 

.3 

.300001 

Inclination 

0® 

0® 

Peri focus Passage 

2,000sec 

2,000sec 

Argument of Perigee 

0° 

0® 

Ascending Node 

0® 

0® . 


Final Integration Time = 10,000 sec 
Again the subtraction of r^ - "Fg can result in a loss of up to 5 
significant digits. The primary differences between this and the previous 
Pair Is the larger eccentricities. Thus enabling another look at the 
effects of a greater eccentricity on the orbits. Two more pairs were 
considered, both pairs were very non-similar, the first one is presented 
below. 

Pair G 

In this case, two near earth orbits which are very dissimilar are 
modeled. The semi -major axis of satellite 2 is chosen such that the ratio 
of the periods of the two satellite is 1/3 resulting in a small separation 
every 3 revolutions of satellite 1. This provides another chance to study 
the propagation of the errors incurred in the first subtraction. 



Figure 7a 

X v&» y for Fair 6 
(KILOMETERS) 


1.00 80.00 120.00 

TIME *10*^ Seconds 

Figure 7b 

TIME HISTORY OF REUTIVE DISTANCE 
FOR PAIR G 
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Orbital Parameters 

Satellite 1 


Semi -major axis 7378.165km 

Eccentricity *3 

Inclination 0® 

Peri focus Passage 2,000sec 

Argument of Peri focus 0° 

Ascending Node 0® 


Final Integration Time = 10.,000sec- 


Satellite 2 
15,347.2km 

.6634883367 

0 ® 

2,000sec 

0 ® 

0 ® 


The minimum separation occurs at perigee and approximately 5 signifi- 
cant digits can-be lost in the subtraction. The pair model resembles Pair D, 
except that the eccentricities are greater for both satellites. The orbit 
shapes and relative vector time history are shown in Figures 7c and 7b. 

The last pair model is again a nonsimilar pair, but they are more similar 
than G or D. 


Pair H 

These again model two near earth orbits, but which are more similar 
than Pair G. The ratio of the periods is now 1/2; and there is one close 
approach every two revolutions of the first satellite. The orbit shapes 
and relative position vector time history is shown in Figures 8a and 8b 
respectively. 
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Semi -major axis 
Eccentricity 
Inclination 
Perifocus .passage 
Argument of Perigee 
Ascending node 


Orbit Parameter 

Satellite 1 
7378.165km 

.3 

0 ® 

2,000sec 

0 ® 

0 “ 


Satellite 2 
10,917.672km 

.526949 

0 ° 

2,000sec 

0 ® 

O'* 


Final Integration Time = lO.OOOsec 
The subtraction of - Fg can again result in a loss of up to 5 
significant digits. 

Once the satelli.te pairs to be modelled have been decided upon, the 
next step is to choose the comparison method and criteria for comparison. 


3.2 Testing Procedure 

The first step in evaluating the relative merits of the various 
formulations is to develop viable comparison criteria. Since, in theory, 
it should be possible in any integration to achieve any reasonable accuracy, 
given a small enough step size, the criterion chosen was efficiency. The 
efficiency of' a method can be based either on the total amount of computing 
time required to achieve a certain accuracy, or upon the number of function 
evaluations (calls to the derivative subroutine) required to achieve a certain 
accuracy. The total integtation time required by each method is obviously 
going to be dependent on many factors such as: 

1) differences in the coding of the method 

2) amount of information calculated at each step for each method, 

3) relative difficulties involved in obtaining certain types of 
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information from each method. 

One instance of the third problem is that if acceleration, velocity, 
and position information for each of the satellites and the relative vector 
are desired in the coordinate system where time is the independent variable 
at each output point, then transformations from s to t are required for 
each of the time smoothed equation sets, while the information is readily 
available for the Cowell formulations. Since the type and frequence 
of information required is very dependent upon needs of the particular 
user, the number of function evaluations required to achieve a disired 
accuracy was chosen as an efficiency criterion. For the particular force 
model chosen (two body only) the function evaluations are relatively 
inexpensive; but for more complicated force models they can become the 
dominant user of computing time. It is felt that trends apparent in this 
study will be reflected in those with more complicated force models. 

For any comparison, a method of determining relative accuracy is a 
necessary consideration. Three approaches for determining relative accuracy 
were examined. 

In the first approach, perturbed motion of the two satellites was 
integrated backwards in time from some reference state to some initial state 
using an integrator which had control over truncation error. The satellites 
were then integrated forward to the reference state, and a bound on the 
integration error was obtained by comparing the final solution to the 
reference state. 

This approach had two major drawbacks. First, a considerable amount 
of setup time was required for each problem, to determine the I'true" or 
reference solution. Second this method would only examine the effects of 
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a particular perturbation history. It was felt that it would be more appro- 
priate to examine the basic behavior of the equation sets,’ so this procedure 
was abandoned. 

Another possibility is to transform the problem into one in which the 
initial conditions could be recaptured at a- given time. The integrated 
solution would then be compared with the initial conditions. This recapturing 
of initial conditions severly restricts the type of problems which could be 
studied, and was discarded in favor of the following method. 

A model was developed with strict two body motion. Such a model allowed 

use of the known integrals of the two body equations. Given the initial 

states for each satellite, reference state vectors for each satellite could 

then be computed using the known two body solution at any time. These 
reference states were computed in double precision to obtain the true 
(reference) states as accurately as possible. 

A program was written to compare the accuracies and efficiencies of the 
methods tested. The logic flow within the program is diagramed in Figure 9. 
The desired initial conditons for each satellite, either in rectangular 
Cartesian or orbital elements are first read in, along with the desired 
central body, (earth or moon), the unit of length to be used (kilometers, 
earth radii, or lunar radii), and the output desired (Table 1). All required 
transformations necessary to get the initial states in the proper form for 
each of the methods are then carried out. These transformations are described 
in the Appendix. 

After all input and initialization is completed, the actual integration 
is carried out. Because of the required time matching for the third 
equation set, and the complications this entails, the third equation set 


V 
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is treated first, then the others all follow the same pattern, with the 
exception of the time matching required for set 3. Equation set 3 is first 
integrated to the first desired output point. This is done by integration 
of the first satellite until the time has passed the first output time, 
then the second satellite is integrated until t2 passes t-j , the integration 
is then iterated until the two times are within 10"^ second of one another. 

It was found that a more restrictive match than this often resulted in a 
stepsize in the s system which was so small as to be numerically zero. 

The time state was then calculated and all transformations required to 
compute desired output information was calculated. The output information 
was then computed and stored. 

At this point a note on the method chosen to acheive time matching for 
equation set 3 is- in order. The process of iterating on stepsize requires 
a great number of added function evaluations, and the efficiency curves for 
the method make the formulation look totally unfeasible. Other methods may 
be devised for accomplishing the necessary time matching, more economically, 
but at this stage, the iteration method seems to be a fair approach. The 
added expense of interpolation would not be visible in the comparison, and 
errors incurred in an interpolation process would be difficult to isolate, 
from the errors inherent in the integration. 

After the required time matching, has been completed, the system 
of equations is integrated to the next output point. The process is re- 
peated until the final time has been reached. 

The integration method employed was a Runge-Kutta (7)8 with Fehlberg's 
coefficients. This code employs variable step size and attempts to keep 
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the relative local truncation error at each step less than the desired 

_3 

tolerance. Each case was run at various tolerances between 10 and 
-8 

10 in order to generate the data for the efficency plots. 

The rest of the equation sets are then processed in the same manner 
sequentially. Although it would be more efficient to integrate the equations 
simultaneously and make use of some of the intermediate information which is 
common to two or more equation sets, it was felt that totally isolating 
each equation set from the others would prevent any possible interaction 
by the equation sets. One such interaction is that since a variable step 
Runge-Kutta method was employed to carry out the integration; the smallest 
step size required for any of the equations would be used for- all equation 
sets. 

After the integration of. all eight equation sets has been completed, 
the desired tables and plots are produced. This included the description 
of the comparison methods employed. The next chapter will present the 
results of these tests. 
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OUTPUT INFORMATION 
WHICH CAN BE SELECTED. 

TABULAR DATA 

1. CP Time required by each method 

2. Function evaluations for each method 


PLOTTED DATA 

(Line Printed, Cal comp, or Both) 


Magnitude of Vector 
vs. Time 

Available For • 


Satellite 

1 

Satel 1 i te 
2 

Relative 

Vector 

Position 

X 

X 

X 

Velocity 

X 

X 

X 

Acceleration 

X 

X 

X 

Angular Momenta 

X 

X 


Error in Position 

X 

X 

X 

Error in Velocity 

X 

X 

X 

Error in Acceleration 

X 

X 

X 

Error in Angular Momentum 

X 

X 


X vs. y 

X 

X 


y vs. z 

X 

X 

■ 

Function Evaluations vs. Error 



X 


TABLE 1 












LOGIC FLOW WITHIN COMPARISON ROUTINE 



.Figure 9' 






CHAPTER IV 


Results of Numerical Studies 

Each of the eight satellites pair models described in Chapter III 
were integrated using the eight equation of motion formulations. The- 
information from these runs was first examined in order to determine the 
most efficient choice of n, the time smoothing exponent for eqqation 
sets incorporating time smoothing. Efficiency plots were made for each of 
the time smoothed equations, the normalized error being plotted against 
the number of functions evaluations for various values of n. The most 
efficient value of n is then the one which lies closest to the origin, 

(the one which yields the best accuracy with the fewest function evaluations). 
To determine the effects of the Nacozy-Szebehely method, methods involving 
it were plotted against their counterparts separately. Thus, Equation sets 
4 and 7 are plotted separately, as are 5 and 8, and 2 and 6. Efficiency 
information is presented for the independently time smoothed equation set 
3 in the first pair model only. Due to the large number of function 
evaluations required for time matching, it was not competat-ive with the 
other methods from an efficiency standpoint. It was found that in- 
clusion of set 3 on the comparison graph caused the scale to become so large 
as to make the efficiency information on the other 'equation sets difficult 
to discern. Thus, set 3 data was omitted for all pairs except pair A. In 
all of the plots, all errors plotted are the normalized norms of the 

errors at the final integration time, i.e., error in 1 1/2 

- _ B^true“^cal cul ated^^‘*’-^^true~'^ca1 cul ated ■ ^^^true’^cal cul ated^ J 

I '"true I 

In the efficiency plots, each of the lines is composed of three data 
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points, the lines connecting them are just to show trends, and do. not denote 
imbedded data points. 

The results will be broken down into eight sections, one for each 
satellite pair. The first portion of each section will deal with a 
determination of the best choice for the smoothing exponent in the time 
smoothing equation sets. Once this has been found, the Nacozy-Szebehely 
version of the relative motion equations will be compared with their’ standard 
counterparts. Finally, the most efficient form of each equation set will 

be presented on the final graph of each section. 

As this comparison is based only upon the errors at the final time, 

it is necessary to determine how the equation sets behave as a function of 
time for each Pair tested. A plot depicting the time history of the normal- 
ized error in distance of the various equation sets is included. for each 
pair. 

4.1 Satellite Pair A 

From Figures 10.1 a and b, the appropriate values of n for equation 
set 4, and its Nacozy-Szebehley counterpart equations, set 7, are seen to 
be 1.0 and 1.5 respectively. With these values a comparison of the two 
equation sets (Figure 10.1c) show that the N-S formulation does help the 
efficiency over a large region. Figures 10. Id, and lO.le show that the 
value n=l is most efficient for both sets 5 and 8. The N-S method appears 
to make equation set 8 slightly more efficient than set 5 at 
high accuracies. Figure lO.lg shows. that set 6, (the N-S version of 
set 2) is again slightly more efficent. In general for this pair, the N-S 
method does not appear to make as large difference in the efficiencies of 
the methods. This is probably due to the fact that the smallest separations 
which the two satellites experience is very large, and the number of signif- 
icant digits lost in the subtraction of rj - Fg is small. Looking at 
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FUNCTIONS EVALUATIONS 


Figure 10, Ig 

Equation Set 2 vs. Equation Set 6 
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TIME HISTORY OF ERROR PAIR A 
Figure 10. li 
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Figure lO.lh, it can be seen that equation set 3,' the independently smoothed 
equation set was very inefficient as previously mentioned. Equation .set 6 

appears appreciably more efficient than the standard Cowell equations, 
especially at higher accuracies, and the time smoothed equation set are 
even more economical for highfep accuracids. 

If the time histories of the equation sets (Figure 10. li) are now 
compared, it can be seen that the various equation sets appear to hold 
their positions relative to one another for the time displayed. From this 
observation -one might be justified in the assumption that the efficiency 
comparisons are valid over a reasonable neighborhood of time around the 
point at which the comparison was actually made. 

4.2 Satellite Pair B 

The satellites modelled in this pair have much smaller separations at 
apogee and perigee and the orbits are more similar than those of Pair A. 
Figure 10.2a shows the value of n to have a more pronounced effect upon 
the efficiencies of equation set 4 than was observed in satellite Pair A 
(Figure 10.1a). The most efficient choice of n for this equation set is 
not clear cut and apparently depends upon the accuracies required, but for 
most cases the choice n=1.5 appears optimum. Figure 10.2b is much more 
definitive in showing n=l to be most efficient for equation set 7. However 
comparison with Figure 10.1b shows the equation set to be more drastically 
affected by the choice of n for this satellite pa1>' than the previous one. 
The effect of the Nacozy-Szebehely method upon equation set 4 is much more 
pronounced for this pair, as can be seen from Figure 10.2c. An accuracy 
of 10”^ can be achieved by Equation set 7 with almost 13% fewer function 
evaluations than with its standard subtraction counterpart. Equation set 4. 
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Figure 10. 2h 

Comparison of Equation for Satellite Pair B 
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TIME HISTORY OF ERROR PAIR B 
Figure 70.21 



52 


Comparison of Figures 10. 2d and 10. Id show that this pair's closer 
approach did not seem to affect the choice of n for equation set 5, but 
merely moved the points to the left and up. The same conclusion can be 
reached for equation set 8 by observing the differences in Figures 10. 2e 
and 10. le. The N-S technique does not appear to have a large effect upon 
equation set 5, as can be seen from Figure 10.2f, but has a better effect 
than it did on Pair A (Figure 10. If) resulting in a savings of almost 6 % 
in function evaluations for most of the graph. Figure 10. 2g shows set 
6, (the N-S version-'df equation set 2) to be consistently more efficient .than 
set 2, although this increase in efficiency is limited to about 6%. 

Comparison of this to the same graph for Pair A (Figure 10. Ig) shows that the 
N-S method definitely helps set 2 for Pair B to a greater extent than 
Pair A, 

The comparison of the best version of each of the equation sets on 
Figure 10, 2h shows that a dramatic increase in efficiency can be realized by 
Equation set 7. The same accuracy as the next best set can be obtained with 
22 % fewer function evaluations. The other equation sets seem pretty much 
clustered with one another. For Pair B, then, the N-S technique improves 
the efficiencies of relative equations and the effect of n on sets 4 and 7 
is. more-pronounced. than for Rain A. 

A glance at Figure 10.21- shows that once again the equations sets all 
appear to vary in the same fashion, with time. This observation seems to 
support the belief that the efficiency comparisons could be accurately 
extrapolated to the region around the chosen comparison time. 

The reader should recall that the primary difference between Pair a 
and Pair B is.-thatf-the orbits in Pair B are much more similar and have a 



53 








54 







NORMALIZED ERROR IN RELATIVE DISTANCE 


55 


CO 



Figure 10.3g 

Equation Set 2 vs. Equation Set 6 


FUNCTION EVALUATIONS 
Figure 10. 3h 

Comparison of Equations for Satellite Pair C 
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smaller separation-. The orbits in the next pair,’ Pair C, are also very similar 
and have small separation distances. Pair C is a high altitude Pair with 
both satellites having very small accentri cities. 

i 

4.3 Satellite Pair C 

The efficiency plots show rather erratic behavior for the methods in 
high accuracy regimes for this satellite pair. It is felt that this is 
due to the effects of entering the region of round-off error. 

Figure 10.3a, when compared with Figure 10.2a, shows the effect of 
n to be much less on the efficiency of equation set 4, but the best 
choice for n appears to be 2.0. In Figure 10.3b, the same type of behavior 
is exhibited for set 7, with n = 2.0 appearing to be the most efficient 
value. The comparison of equation sets 4 and 7 in Figure 10.3c shows the 
N-S method to be very advantageous for this equation set. Note that a 
savings of up to 25% in function evaluations can be realized over most of 
the range of accuracies. The determination of the best n for sets 5 and 8 
from- Figures 10.3d and 10. 3e, yields n = 1 to be the best but it should 
be noted that these equation sets are not as clearly affected by the value 
of n as were equation sets 4 and 7. 

Figures 10. 3f and 10. 3g show little benefit from the application of the 
N-S technique to sets 5 and 2. This is a major change from the previous 
pair. The comparison of the most efficient version of each equation set, 
presented in Figure 10. 3h, does not show any one set to be most efficient 
for all accuracy ranges of accuracies, but shows set 7 to be the best over 

I 

the largest range of accuracies with the other sets clustered near it. 
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FUNCTION EVALUATIONS 


Figure 10. 4g 

Equation Set 2 vs. Equation Set 6 
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Figure 10.4h 

Comparison of Equation for Satellite Pair D 
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Iri' stiiwary, for >P4ir ' C', which primarily differed dniy in eccentricity 
from the Pair B-, the N-S method does not seem to help as much, and n ' 
not affect the time smoothed equations as drastically as before, although 

for sets 5 and 8 the exponent n=il now appears best. 

Figure 10.31 shows the equation sets to vary little relative to 

one another as the time increases. The calculated efficiencies at the final 
time the, would appear to be representative for this pair. 

The next Pair models two very dissimilar orbits, which differ even more 
than Pair A. 

4.4 Satellite Pair D 

Figures 10.4a, b, and c are very similar to their counterparts for 
Pair A, in that the value of n affects the efficiencies of equation sets 
4 and 7 to a lesser degree than in pair B. Again, for set 4, the most 
efficient choice of n appears to oscillate between 1 and 1.5, with 

n = 1.5 being better over a larger region than for n = 1.0. 

The value of n for set 7 also has small effect again, but this time 

choice of the most efficient value is not clear cut. However, n = 2 
seems to be a good compromise, due to its high efficiency at greater 
accuracies. Figure 10.4c shows the N-S implementation in set 4 does not 
improve matters but is very close to set 7, as in Pair A. 

Figures 10. 4d, e, and f show sets 5 and 8 to be relatively insensitive 
to the value of n, but the N-S implementation seems to have helped slightly. 
Figure 10. 4g displays the change in efficiencies of set 2 with and without 
the N-S technique, and again it doesn't affect the accuracy or efficiency 
much. 

Figure 10. 4h compares the most efficient foms of each of the sets, 
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and set 1 appears much better than the others. In summary, for these very 
dissimilar orbits, the N-S method does not appear in general to be as 
effective as the striaght subtraction, (standard Cowell). In an attempt to 
understand why set 1 is so much more efficient, it may be advisable to look 
at the relative distance vs. time plot for this pair (Figure 4b). As the 
satellites approach the point of close separation, the subtractions needed 
by the methods which require integration of the relative vector can cause 
a very abrupt loss in precision. This presents problems for the other 
methods and they must reduce step size 

problem and can continue at a larger step size. Another important aspect 
of this pair can be seen from Figure 4b. Due to the final time chosen, 
the satellites have a large separation when the relative vector is calculated 
for equation set 1 , and a minimal loss of significant digits occurs in the 
subtraction of ~ ^2* standpoint, the fact that the time 

smoothed equations which integrate the relative vector competed as well as 
they did is quite impressive. 

Figure 10. 4i shows the equation sets to be reasonably well behaved, 
however there is quite a bit of range in the relative efficiencies of the 
methods at different time points, and no extrapolation to times around the 

comparison points can safely be made. 

The effect of the altitudes on the equation will be studied next by 

comparing the results of Pairs C and E which except for their altitudes 

are almost indentical. 

4.5 Satellite Pair E 

Figures 10.5a, b, and c are almost identical with those from Pair C 
(for Equation Sets 4 and 7) in shape. However, when Figure 10.5c is 
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compared to 10.3c it can be seen that the relative effectiveness- of the 
N-S method on set 4 has decreased with decrease in altitude. 

Figure 10. 5d is very similar to its equivalent in Pair C, while 
for equation set 8, Figure 10.5e, the effect of n on the methods is 
much more pronounced. In general, from Figure 10. 5f the efficiency of 
equation set 5 is greatly improved by the N~S technique, while for Pair" 

C there was not much effect. 

Equation set 2 also is greatly aided by the N-S technique, as can be 
seen from Figure 10. 5g. The comparison of the best form of each equation set. 
Figure 10.5h, shows set 8 to be the most efficient over most of the range. 

The major effect of the decreased altitude between Pairs C and E 
was to increase the amount the N-S technique helped equation set 5, while it 

did not seem to increase the efficiency of equation set 4. 

In an attempt to see if the efficiency comparisons made at this 

particular time can be extrapolated for other time periods within the 
the integration, a study of Figure 10. 5i is in order. The equation sets 
oscillate drastically in the early part of the integration, then become 
more smooth near the comparison point. Even near the comparison point 
there is some movement of the sets relative to one another, and it is 
apparent that the relative efficiency comparisons are only valid at the 
point compared. 

A comparison of Pairs E and F will help to show the effect of changing 
the eccentricities on the equation sets. 
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4.6 Pair F 

A comparison of Figures 10.6a and b with 10.5a and b tend to show that 
an increase in eccentricity has only a small effect on the relative 
efficiencies of equation sets 5 and 7 at various values of n. In contrast, 
Figure 10.6c shows a dramatic increase in efficiency realized from the 
Nacozy-Szebehely method applied to equation set 4. However, this efficiency 
increase is not as drastic as it was for Pair E. Figures 10. 6d and e 
vs. 10. 5d and e again show little change in the effects of the value of n 
for Pairs E and F, sets 5 and 8. Again, the increase in efficiency pro- 
vided by the N-S technique has increased for set 5 with the increased 
eccentricity. 

Figure 10. 6g shows a dramatic gain in efficiency when- the N-S method 
is applied to set 2, (set 6), especially at high accuracy regions. 

Figure 10. 6h shows the two time smoothed equations compared to their 
performance in Figure 10. 5h h ave improved a great deal relative to the 
nonsmoothed equations. However, set 6 has gotten a little less efficient 

although it is still better than set 1 over most of the range. 

Figure 10. 6i displays the time history of the error norm used for 

each equation set. The equation sets are all well behaved, and do not 
cross one another, this tends to show that the exact time chosen for 
comparison would not critically affect the conclusions drawn. 


4.7 Pair G 

Pair G, like Pair D, approaches the point of close separation abruptly. 
The N-S method does not increase the efficiency of set 4 (Figure 10.7c), 
while it is competative for set 5 vs. set 8, (Figure 10. 7f) and offers 
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a decided advantage for equation sets 2 and 6, (Figure 10. 7g). The 
integration of the relative vector (Figure 10. 7h) does not help the efficiency 
of the Cowell methods, while the time smoothing does seem to offer some 
advantage, especially for set 5. Again, a look at the magnitude of the 
relative position at the final time Figure 7b shows that Equation set 1 was 
not hampered by having to do a subtraction of two nearly equal numbers 
at the final term, and yet equation set 5 which must calculate relative 
accelerations by subtracting two nearly equal numbers (when satellites are 
at perigee) is still more efficient; This would tend to substantiate the 
feeling that the smoothing process in itself helps the efficiency of the 
methods . 

The importance of the time chosen for the comparison alluded to above 
can be easily seen from Figure 10. 7i. The relative accuracies of the 
equation sets varies widely, and a comparison at another time might show 
quite different results. 

4.8 Satellite Pair H 

From Figure 8b it can be seen that this pair model also approaches 
the point of minimum separation abruptly, but not as sharply as Pair G. 

The effects of this on the use of the N-S method are to make them more 
efficient for sets 4 vs. 7 (Figure 10.8c) and very competitive for sets 
5 vs. 8 (Figure 10. 8f). On the Modified Cowell sets (set 2 and 6), the N-S 
equations become more competative as the accuracy requested increases. But 
it appears that set 2 is more efficient at low accuracies. Figure 10, 8h 
shows the integration of the relative vector in the Cowell methods to be 
less efficient, while the use of time smoothing appears to increase the 
efficiency, especially for set 8. 
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The reader's attention is now directed to Figure 10. 8i, the time 
history plot for Pair H . It should be observed that once again the 
relative accuracies of the methods vary drastically with time, and for 
this reason any conclusion drawn from the efficiencies of the sets at the 
comparison may not be valid for other times. 

4.9 Data Comparisons 

Since the needs of each individual user of this data will be dif- 
ferent, it would be impossible to anticipate all required comparisons, 
and make them in this report. However, since all the data are presented, 
the user can make further comparisons as needed for any purpose. Selected 
comparisions will now be presented which illustrate the information 
available in the data. 

The next step is to determine how Equation set 5 (time smoothing 
applied to one satellite and the relative vector based upon the magnitude 
of both satellites position vector) reacts to different situations. For 
the case where the closeness of approach is varied (Pair A vs. Pair B) 
it can be seen from Figures 10. Id and 10. 2d that the efficiency of set 5 
is very sensitive to changes in n for both satellite pairs, and n = 1.0 
is most efficient by far for both. A comparison of Figures 10. If and 10. 2f 
shows how set 5 relates to set 8, which is plotted in Figures lO.lh and 
10.2h. This comparison shows that while for Pair A, set 5 is very competi- 
tive with set 4, and more efficient than either the Standard or Modified 
Cowell procedures, it is clearly less efficient than set 4 for Pair B and 
has moved much closer to the standard Cowell procedure. From this it would 
appear that the introduction of the magnitude of the second satellite into 
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the smoothing equation has decreased the gains realized by the integration 
of the relative vector. 

The reaction of set 5 to a change in. semi -major axis is next examined 
through a comparison of Pairs C and E. Figures 10.3d and 10. 5d show set 
5 to be sensitive to changes in n, but the most efficient choice of n 
in both' cases is still n = 1.0. When the plot for set 5 is compared to 
set 8 in Figure lO.Bf, and its location extrapolated relative to set 8 
in Figure 10. 5h, a comparison with Figure 10. 3h shows that the change in 
semi -major axis has not had a noticeable effect upon set 5, relative to 
sets 1, 2, or 4. The effect of the change in eccentricities between Pairs 
E and F is first noticed in Figures 10. 5d and 10. 6d. Set 5 is still 
sensitive to changes in n, but n = 1 is most efficient for both pairs. 
Again, in order to compare set 5 to sets 1,2, and 4 it is first necessary 
to look at Figures 10. 5f and 10. 6f to note its position relative to set 8. 
Once this is done, it can be seen that although in both cases the efficiency 
is better than that attained by set 2, it is again hampered in the 
high eccentricity case {as was the other method involving integration of 
the relative vector). Set 5 and set 4 are very close to one another although 
it appears that in the case of the low eccentricity (and therefore less 
abrupt approach) set 5 was more efficient while for Pair F, set 4 was 
more .efficient. 

To further study the effect of the abrupt approach on Set 5, Pairs D, 

6, and H are now compared. Figures 10. 4h, 10. 7h andl0.8h all exhibit 
the same trend. Over the entire range, set 5 is more efficient than sets 
4 and 2, and for Pairs 6 and H, set Sis much better than the Standard 
Cowell, but seems to be approaching it as the accuracies get lower. For 
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Pair D, set 5 is not quite as efficient as the standard Cowell. The 
overall conclusions are in general that set 5 has better efficiency than 
the other methods involving the integration of the relative vector. However, 
it suffers from the limitations imposed upon the methods by the abruptness 
of approach of the two satellites. 

In order to observe how the integration of the relative vector is 
affected. by the variation in the closeness of approach between Pairs A and 
B it should first be noted from Figures 10. Ig and 10. 2g that the N-S method 
does not affect the set 2 equations drastically. A comparison of lO.lh and 
10. 2h shows, however, that the Modified Cowell method with the N-S implemen- 
tation (set 6) is more efficient than the Standard Cowell formulation 
(set 1) in both cases. Note that for PairB, with its much closer approach, 
the relative gain in efficiency is far greater. 

The way that a change in the semi -major axis of the two satellites 

affects the relative effi cienc y of the Modified Cowell equation can be 

investigated by comparing the results of Pairs C and E. Figures 10. 3g and 

lO.Bg show a larger change in efficiency for set 2 with the implementation 

of the N-S technique. Figure 10. 3h and 10. 5h both show the Modified Cowell 

with N-S technique to be competitive with the Standard Cowell method for 
-6 

low accuracy (10“ and above) regions, while in high accuracy range 
(below 10“^) the Modified Cowell procedure is much more efficient. On both 
pairs, almost one order of magnitude gain in accuracy can oe achieved with 
the same number of function evaluations. From this comparison, it appears 
that the size of the orbit does not significantly affect the relative 
efficiencies of the Standard Cowell and Modified Cowell formulations. 
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The next effect studied is how changing the eccentricities of both 
satellites affects the efficiency of the Modified Cowell equations relative 
to the Standard Cowell set. Pairs E and F are almost identical except 
that the eccentricities of both satellites are increased, while this will 
change the size of the orbits also) from above we have seen that this 
effect is slight). The first thing to note -is that in Figures 10. 5g and 
10.6g, the application of the N-S method has increased the efficiency of 
set 6 by an appreciable amount. Now, from comparision Figures 10. 5h and 
10.6h, it can be seen that for both Pairs E and F, the Modified Cowell 
are again competitive with the Standard Cowell for low accuracies. For- 
high accuracies, in Pair E the modified Cowell set is much more efficient 
than the Standard Cowell, while for Pair F the reverse is true. This 
comparison would tend to substantiate the observation that as the orbits 

f 

become more eccentric, the advantages gained by integration of the relative 
vector become smaller. A comparison of the time history of the magnitudes 
of the relative positions of the two satellite pairs (Figures 5b and 6b), 
suggests a possible reason for this. From the figures it can be seen that 
for the more eccentric pair. Pair F, the satellites approach one another 
more abruptly at perigee (t = 22 seconds and t = 85 seconds) than the 
near circular orbits (Pair E). The separation pf the two satellites for 
these two pairs is approximately the same, while the eccentric satellites 
are nearer the central body. 

For the Standard Cowell Set, this abrupt approach is no problem as the 
integrator will just integrate the orbits independently and not know of the 
change. Figures 5b and 6b also show that at the final time where the errors 
areccalculated, the relative distance is quite large, so the calculation of 
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F = Tg - results in the least loss in significance for set 1, and is 
a very optimistic situation for the set. 

To -further observe the effects of this type of abrupt approach, the 
reader's attention is directed to Pairs D, G, and H. From Figure 10. 4h it 
can be seen that, as expected, the integration of the Modified Cowell 
set is not as the Standard Cowell set. From an examination of Figure 4b, 
the abruptness of the approach can again be seen. Figure 10. 7g shows 
that the implementation of the N-S technique has greatly improved the 
efficiency of the Modified Cowell set, but even so, from Figure 10.7h, set 
6 is much less efficient than the Standard Cowell set. Again, the time 
history of relative position for Pair E, (Figure 6b) shows a very abrupt 
approach. Finally, Figure 10. 8h shows this same trend; the integration 
of the relative vector in set 2 has drastically reduced the efficiency of 
the Cowell equations. Figure 8b shows the abrupt approach which hampers 
the relative equations. 

The integration of the relative vector in general increases the 
efficiencies of the Cowell type of equations when there are no abrupt 
approaches, and the benefits of this Modified set appear to fall off as the 
approaches ‘become more abrupt. The trade-off point for efficiency between 
the Modified Cowell and. the Standard Cowell sets is dependent upon 
how rapidly the relative vector changes near the close approaches. 

The behavior of equation set 4 (time smoothing applied to one satellite 
and the relative vector based upon the magnitude of the satellite position 
vector) will now be discussed. The first effect studied is how the accuracy 
of the equation is affected by a change in the closeness of approach 
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between Pairs A and B. As pointed out previously, the value of n has a 
more dramatic impact upon Pair B for set 4 than for Pair A. For Pair A, 
the largest change in efficiency is between n=l and n=2 and is approx- 
imately 30%. The other major difference between these curves is that for 
Pair A, n=l is most efficient for most accuracies, while for Pair B 
n=1.5 is better over most of the range. From Figure 10.1c, it can be seen 
that the inclusion of the N-S technique does not improve accuracy for Pair A 
a great deal. From Figure 10. 1h it can be seen that set 7 is more 
efficient for Pair A than the Standard Cowell method. From a comparison 
of Figures 10.1c and h, it can be seen that for Pair A, set 4 is more 
efficient than the Standard Cowell set as well as the Modified Cowell set 
at high accuracies, and .competitive with them for lower accuracies. 

If the graph of Equation set 4 vs. set 7 (Figure 10.2c) is now com- 
pared with Figure 10.2h, it can be seen that set 4 for Pair B is much more 
efficient than the Standard Cowell set, and the Modified Cowell set. The 
relative improvement in set 4 between Pairs A and B is dramatic. This 
improvement is even greater than that realized by integration of the 
relative vector rather than the Standard Cowell set. Thus, the time 
smoothing has apparently greatly improved the efficiency of the integration, 
and this is not due to the integration of the relative vector alone. 

To see how the semi -major axis size affects the efficiencies of set 4j 
examination of Pairs C and E is useful. First', from Figures 10.3a and 
10.4a it is clear that this change in size has not changed the relative 
effects of n on the set. If the reader looks at Figure 10.3c to see 
how set 4 looks relative to set 7, then looks at Figure 10. 3h to approximately 
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locate set 4 relative to set 7, and then follows the same procedure for 
Figure 10.5c and h, it can be seen that the efficiency of set 4 relative to 
set 1 for Pair E has not been appreciably effected. Thus, the efficiency 
gain realized by the time smoothing does not appear to be significantly 

effected by the orbit size. 

/ 

The effect of the eccentricity on set 4 can be observed in its 
efficiency relative to the Standard Cowell equations for Pairs E and F. 
Figures 10.3a and 10.6a (difference is eccentricity) show that set 4 is 
not extremely sensitive to the value of n for either Pair. Again, in 
order to compare the relative efficiency of set 4 relative to the Cowell 
and Modified Cowell equations it is first necessary to observe the curves 
of Set 4 relative to set 7 for Pairs E -and F. 

Next, compare set 7 to sets 1 and 2 in Figures 10.5c and 10.6c. 

Once this is done, it can be seen that for both pairs, the time smoothing 
has increased the efficiency beyond that for set 1 for both pairs in the 
low accuracy region. However, for Pair E in high accuracy regions, time 
smoothing results in a more efficient procedure than set 1 while it is 
almost the same as set 1 for Pair F. It should =be observed that this is 
the same trend as was observed for the Modified Cowell set, and indeed 
•this set also integrates the relative vector. 

If set 4 is compared to set 2, it appears that set 4 is always more 
efficient than set 2, due to the time smoothing, while set 4 reacts in the 

same manner as set 2 in the case of an abrupt approach. To substantiate 
this, a comparison of Pairs D, G, and H is in order. Comparison of Figures 
10. 4h, 10. 7h, and 10. 8h tend to support this hypothesis. From the figures 
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the smoothing does tend to increase the efficiency relative to set 2, but 
the fact that the relative vector is being integrated hurts the efficiency 
of the method relative to set 1. 

In summary, it appears that the time smoothing enhances the efficiency 
of the integration of set 4 a great deal. However, the integration of the 
relative vector introduces the same form of problem for this set (set 4 as was 
noted for set 2, the Modified Cowell set, when abrupt approaches are 
integrated. Again it is necessary to remember that due to the final time, 

I 

the efficiencies given for set 1 are very optimistic, and would likely 
be much lower if output were desired at a point of minimum separation. 

The Nacozy-Szebehely method applied to sets 2, 4 and 5 (Sets 6, 7, and 
8 respectively) had very consistent effects. If the differences between 
each of these sets is looked at for Pairs A and B, it can be seen that 
the N-S technique increajsed the efficiencies of each of the methods for 
higher accuracies for Pair A, and increased their efficiencies even more 
for the closer approach in Pair B. 

The semi -major axis difference between Pairs C and E has an interesting 
effect upon the efficiencies of sets 2, 4, and 5 and their N-S counterparts. 
For Pair C, there is little increase in efficiency for the N-S versions 
at high accuracies, while for Pair E there is a dramatic increase for 
each of the methods at high accuracies. For the comparison of Pairs E 
and F, the implementation of the N-S technique behaves similarly for both 
However, it is even more efficient for the moderately eccentric case than 
for the near circular case. 
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The Increases in efficiency due to the N-S procedure in the above 
cases are quite large, in many cases allowing an accuracy increase of 
almost two orders of magnitude for the same number of function evaluations. 
A comparison of Pairs D, G,*and H shows a little more erratic behavior in 
the advantages of the N-S method. For Pairs D and G, the method has little 
effect on sets 4 and 8 (Figures 10.4c and 10.7c) but tends to decrease the 
efficiency, while for Pair H, the efficiency is increased slightly. 

Equation sets 5 and 8, on the other hand, are slightly aided by the N-S 
technique in Pair B, while it is slightly hampered for Pairs G and H. For 
Pairs D, set 2 is only slightly effected by the N-S method, while for 
Pair 6 there is a dramatic increase in its effi,ciency, and finally for 
Pair H there is only a small decrease in efficijency. The N-S technique 
appears to help the methods when there is a sustained close approach, but 
when the approach is rapid and short, the method will in general make 
little difference. 

While it is hoped that this chapter has presented most of the major 
comparisons, the graphs included will allow the reader to make any other 
comparisons which he may desire. The conclusions reached will now be 
grouped in a more compact form in the next chapter. 



CHAPTER V 

Conclusions and Recommendations 

It is concluded that for the following cases the Standard Cowell set 
will either be most efficient, or very competative with the other sets: 

1) The desired accuracy on the relative vector is low 

-3 

(10 . and above), or 

2) The satellites have an abrupt approach. 

For the time smoothed equations, the first factor which it was 

' 9 

necessary to determine was the fact value of n. As stated by Nacozy , 

the most advantageous choice for n, when studying one vehicles motion is 

dependent upon many factors, such as 

1) location in the two body orbit, 

2) the order of the integration method, 

3) forces included in the perturbation model, and 

4) ' the eccentricity of the two body reference orbit. 

Clearly, the most efficient choice for n for each equation set cannot be 
determined for all cases, but will only be discussed for the cases run 
here. 

The proper choice for n, the time smoothing exponent for equation set 
4 varies, from pair to pair, and it is interesting to note that the 
introduction of the N-S. method to perform the subtraction can alter the 
most efficient choice of n (i.e., set 7). The proper choice of n seems- 
to be very problem dependent, and will probably have to be determined for 
the specific model being used. The value n=2.was most commonly best for 
sets 4 and 7. 
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Equation set 5 was very consistent in that n=1.0 was most efficient. 
This behavior was repeated in its N-S counterpart, set 8. It is interesting 
to note that in the time transformation function, for sets 4 and 7 

f - l/(r2>" (29) 

and the transformation function for set 5 and 8 

f = l/(r^r2)" (30) 

the choice of n=2 for Eq. (29), and n=l for Eq. (30) produce values 
which are very close to one another if 

The N-S method appears to enhance the efficiencies of the methods in 
almost all cases. This is most evident when the satellites are in close 
proximity for a reasonable length of time, i.e., the orbits are similar. 

In conclusion it is recommended that when the orbits are dissimilar, 
that is, they approach the point of close separation very abruptly, the 
standard Cowell equations will probably be either most efficient, or very 
competative with the other methods described here. When the two orbits 
are similar, the modified Cowell equations, incorporating the N-S technique 
will probably offer great advantages at high tolerance. A very important 
point to note here is the fact that the machine used (CDC 6600) carries 
approximately 14 digits, as the close approaches in which 5 digits were 
lost will probably not show the advantages of the N-S technique as clearly 
as on a machine of lesser word length. Another point to be considered is 
that the N-S method requires extra time to do the subtractions. cases 
where the force model is quite simple, the time savings of this method 
might be reduced, but in cases when a more sophisticated force model is used. 
Its advantages will be increased. The time smoothing equations can be 
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very economical to use, but the amount of trajectory-time information 
required by the user and frequency with which this information is required 
can reduce its overall economy as far as computing time is concerned. If 
only position information is required, then no transformations between t 
and s are required after the initialization. If velocity and acceleration 
data is required, however, the transformations might become too expensive 
if the force model is simple. If the force model is more complex, the 
amount of time spent in the transformations may be small compared to the 
savings realized from fewer function evaluations. Again, it must be pointed 
out that these conclusions rest on the basis that the trends in relative 
efficiencies of the methods desired using a two-body model carried over 
when more sophisticated force models are used. Since the two body acceler- 
ations are dominant in most instances, it is believed that while the optimal 
value for n will probably change for each model, the overall trends 

will extrapolate the more complex force models. 

Equation set 3 suffers from the effects of the method used for time 
matching. If a much more efficient method can be developed for determining 
information for both satellites at the same time, such as an interpolation 
procedure, it is felt that this- set might be competative with the other 
time smoothed sets, and for the case of dissimilar orbits might be more 
efficient, (as the Standard Cowell is more efficient than modified Cowell 
for these orbits). 

This paper scratches the surface of a class of problems which 
deserve greater study. Other topics of study which could be investigated 


are: 
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1) An investigation of the effects of the word size of the 
computer used on the efficiencies realized in the implementa- 
tion of the N-S technique. 

2) More efficient means of obtaining information about both 
satellites from equation set 3. 

3) The possibility of choosing the best value for the time smooth- 
ing exponent at each step, analogous to the method proposed 

by Nacozy for Integration of a single body. 



Appendix 


As stated in Chapter Z, the time smoothed equations are all based on 
a transformation of the independent variable of integration, from t to 
a new independent variable of integration s. This transformation usually 
adds a new differential equation to the system, relating t and s. 

The common form of this relationship is 

I 

ds = f dt (1) 

where f is a function of magnitudes of the position vectors of satellite 
1 and/or satellite 2. 

The first step in obtaining the time smoothed equations is to determine 
the relationship between the acceleration of a vector F with respect to 
t, F » and with respect to s, F". The second derivative of J with 
respect to time can be written 



If ~ is replaced by f as indicated by equation (1) and the differentia- 
tion is carried out, the equation can be written as: 

A r^i f + [fi 
dt [ds J ds dt □ 

dt ds dt [ds J • ds dt [ J 
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Letting n denote differentiation of a vector n with respect to 
time, and iv' denote differentiation with respect to S, equation (3) 
becomes 

t = ?" f ^ + ? f (4) 


This equation relates accelerations in the time space to accelerations in 
the time smoothed space. The corresponding relationship between t and s 
can be obtained from equation (1) by a separation of variables and integration. 
The result is 



t = 



ds 

f 


(5) 


In the first time smoothed equation set, each satellite is independently 
smoothed with respect to the magnitude of its position from the central 
body. Letting 

ds^ 1 

^ dt^ ®l*"l*^^ 
and 


ds« 1 

f ss ss 

2 n V. Til 

dtg a-jr^ 1 

it can be noted that f-j and fg 
For this form 


are both of the form’ 


f - 1 

ar 


( 7 ) 
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f = d? 


-1 


“2 


= -(ar*^) na — 


or finally. 


-n ^ _ -n dr ^ 
ar"*’ ds dt-^^n+1 


f = 




dr 

ds 


( 8 ) 


Therefore , 


^1 




-n 


1 




ds. 


1 


a^r]^”2 


f - ~"2 '"'2 




1 


a2.f2 


(9) 


Substitution of r^ and V 2 individually into equation (4), and the 
application of equations (9) to the result yields 


•'l = 


1 


a 2ni+l 
®1 ^1 ‘ 


and 


'‘i = 




dr^ dr-.j 
api^*^!^^ dSi ds^ 


n2 dr2 drg 


^Z'^Z ^ 


ds^ ds« 
S 2 


( 10 ) 


Letting the accelerations in time space be expressed as: 
«. ~u(K|) 


"l = 


3 ^ P'1 (i^l >1''] ) 


and 


>"2 = 


■1 

-u(F2) 


+ P2(r2»r2) 


( 11 ) 


the accelerations of r^ and can be introduced into equations (10). 
Solving the resulting equations for F" to obtain the equations of motion 
for fime smoothed space the equations become 
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„ n, dr, dr, 

y* H s _J_ L 1 

1 r^ ds^ ds.j 


- u a. 


r 7 

n ' >^1 


+ P (r 


,r,) 


.2 2n 
■®1 *^1 


and 



dr2 drg 

r^di^di^ 




( 12 ) 


At this point it is noted that since each satellite has been independ- 
ently smoothed, the relations between t and s are different for each 
satellite. If f^ is substituted into Equation (5), the relationship 
between t and s for the first satellite is found to be 

•^1 “ C* '■1"' ‘‘"l (13) 

no 

While for the second satellite the relationship is 


I ^^2 04) 

^20 

Thus for this equation set, two new first order differential equations have 
been introduced. A major disadvantage of this formulation is that since 
relative position is our primary interest, it will be necessary to match 
t-j and tg at each desired output point. 

Similar logic as was used in developing the Cowell equations for one 
satellite and the relative vector may be used to develop a second set of 
time smoothed equations. For this set, (Eq. set 2) one satellite and 
the relative vector are smoothed, both with respect to the magnitude of 
the position of the second satellite, {can be 1 or 2). The relation 
between t and s is, as before, 

f = 1 

dt 821^2*^2 


(15) 
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The expressions for f and f derived previously and stated in 
» " ♦♦ 

equations (9) are still valid and need only be applied to F and 

F 

c 

it « 

2 


2 . The equation for r 2 is unaltered and is restated below 


rg dSg dS2 ^ ’^2 ^ ^ 


(16) 


The first step in developing the expression for r" is substitution 
into equation (3). 


r = r" f^ + r* f (17) 

2 

If the equations for f and f are substituted into equation (17), and 
F is replaced by 


r = ri - r2 = -V 


1 

V.3 ■ „ 
[r] r, 

the following equation is obtained: 


+ P(r^»r-,) - P(r2,rg) 


( 18 ) 


P(r^.r^) - P(r2,r2) +w 


"2 " 1 ^ 


3 " 3 


1^2 ’"IJ 

- 


(19) 


It should be noted that one direct consequence of this coupling is to 
eliminate the need for two relations between \ and s. The time can be 
directly calculated from equation (5) as 


t = ^ (20) 

•'so ^ 

Now the time matching required for the first set of time smoothing equations 
is not necessary. 

A third set of smoothed equations, suggested by Dr. V. G. Szebehely 
can be obtained by smoothing one satellite and the relative vector with 
respect to the product of the magnitudes of the positions of both satellites. 
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Here f has the form 


" dt 


-n 


( 21 ) 


and f can be developed as follows 


-n-1 


f = -nCr^rg) r^ 


dr« dr,l 

_^ + r — i 

dt ^ dt 


= -n 


-(n+1) f dr, ds dr, ds 

r, ~ + r, — ^ —I 

' ds dr ds dt 


-(n+1) 

= -n (r^rg) • f • (r^rg' + r 2 r^ ) 


= -n (r^r 2 ) 


-(2n+l) 




( 22 ) 


As before, the expressions for f and f can be applied to equation (4) 


for r 2 » resulting in 


1 


2 2n 

(r,r2) 


(r,r,' + r,r ■) 


(r^rj) 


2n+l ^'T2 ’2'1 


(23) 


Solving for r 2 " one obtains 


_ ^ 2n n _ 

r2"=r2(r,r2) + (r,r2' + r^r, ■) 


(24) 


When the acceleration of r 2 with respect to time from equations (11) 
is introduced into equation (24), the smoothed equation for F 2 emerges as 


n - 2n -yr, _ ^ 2n 

'^2"= (r^^2'-*-Vl'^’^2'-("l’^2) -^+P(^2'’"2^("r2^ (25) 

^ 2 r2 


The relative acceleration vector in the transformed system now needs 
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to be developed. Recalling the acceleration from equation (18), 


-r^ r 2 - 

. 4 * I 

r"3 “ ^ 3 

^ I 

\S} ^2 J 

L 


Pl(i^l»ri)-P(r2,r2) 


r = ri - 1^2 ' 


and substituting f, f, and r into equation (17), we arrive at ‘ 


(18) 


-y 


>1 r2 
3 " 3 

+ 1 

h ’^2j 



P(r^r^)-P(r2»r2) 


-2n 


=r“(r,r2) --5-(r,r2-+r2r,') 


r,r2 


Rewriting the equation to isolate r" , produces 


r" = 


1'2 




IL fy^ y» y^ M y* » _ y, y^ ^ + J ! ! ^ ^ ( 26 ) 

r , 2 *^2^ " L >^1 '^2 -1 /. „ ^-2 n 


(r,>-2) 


-2n 


(r,r2)- 


Once again the relationship between time and s may be obtained from 

equation (5) as 

fS-F 

t = J ^ ( rir2 ) ds 

"o 

Again there is no need for the time matching required by the first 


set of smoothed equations. 
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